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Abstract 

The atomic structure and electronic properties of the tip apex can strongly affect the contrast 
of scanning tunneling microscopy (STM) images. This is a critical issue in STM imaging 
given the, to date unsolved, experimental limitations in precise control of the tip apex atomic 
structure. Definition of statistically robust procedures to indirectly obtain information on the 
tip apex structure is highly desirable as it would open up for more rigorous interpretation and 
comparison of STM images from different experiments. To this end, here we introduce a statistical 
correlation analysis method to obtain information on the local geometry and orientation of the 
tip used in STM experiments based on large scale simulations. The key quantity is the relative 
brightness correlation of constant-current topographs between experimental and simulated data. 
This correlation can be analyzed statistically for a large number of modeled tip orientations and 
geometries. Assuming a stable tip during the STM scans and based on the correlation distribution, 
it is possible to determine the tip orientations that are most likely present in an STM experiment, 
and exclude other orientations. This is especially important for substrates such as highly oriented 
pyrolytic graphite (HOPG) since its STM contrast is strongly tip dependent, which makes 
interpretation and comparison of STM images very challenging. We illustrate the applicability 
of our method considering the HOPG surface in combination with tungsten tip models of two 
different apex geometries and 18144 different orientations. We calculate constant-current profiles 
along the (1100) direction of the HOPG(OOOl) surface in the |H| < 1 V bias voltage range, and 
compare them with experimental data. We find that a blunt tip model provides better correlation 
with the experiment for a wider range of tip orientations and bias voltages than a sharp tip model. 
Such a combination of experiments and large scale simulations opens up the way for obtaining 
more detailed information on the structure of the tip apex and more reliable interpretation of 
STM data in the view of local tip geometry effects. 
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I. INTRODUCTION 


The interpretation of scanning tunneling microscopy (STM) images is not straightforward 
due to the effects of the local tip apex geometry, termination and orientation. The reason is 
the convolution of sample and tip electronic states in a given energy window defined by the 
bias voltage, and the fact that in STM experiments the detailed atomic geometry around 
the tip apex is practically unknown and hardly controllable. On the other hand, it is clear 
that the electronic states and their dominating orbital characters involved in the tunneling 
depend very much on the local atomic structure of the tip apex. 


It has been a challenge to obtain information about the relevant properties of the STM 
tip apex for a long time. Herz et al. performed reverse STM imaging experiments to study 


a. 


The 


p, d, and / orbital characters of the tip apex atom above the Si(lll)-(7 x 7) surface 
combination of STM experiments and simulations on well characterized surfaces to obtain 
information on the tip structure and termination was used, e.g., by Chaika et al {2,0]. They 
considered the highly oriented pyrolytic graphite (HOPG) surface in the (0001) crystallo¬ 
graphic orientation in combination with W(OOl) tip models. Rodary et al. studied Cr/W 
tip apex structures by high resolution transmission electron microscopy, and they pointed 
out that the magnetization direction of monocrystalline nanotips cannot be controlled in 


spin-polarized STM |^. Recently, the effect of the tip orbitals on the STM imaging of 
supported molecular structures attracted considerable attention. Gross et al. investigated 
pentacene and naphthalocyanine molecules on NaCl/Cu(lll) surface by CO-functionalized 
tips, and they explained the obtained STM contrast by tunneling through the p-states of 


Q 


the CO molecule ol. Siegert et al. developed a reduced density matrix formalism in com¬ 
bination with Chen’s derivative rule p to describe electron transport in STM junctions for 
molecular quantum dots, and studied the effect of selected tip orbital symmetries on the 


STM images of the hydrogen phthalocyanine molecule on a thin insulating film j^. Lakin 
et al. proposed a method to deconvolute STM images and determine molecular orientations 
of both the sample and the functionalized tip |s|. In their work a C6o-Si(lll)-(7 x 7) surface 
and a Ceo-functionalized tip were chosen. 


Even in seemingly less complicated STM junctions, only a few theoretical works focused 
on the effect of the tip orientation on the STM images. Hagelaar et al. demonstrated that 
a wide range of modeled tip terminations and orientations can reproduce the experimental 
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B- 


images for NO adsorbed on Rh(lll) [9|. This work also showed that the modeling of realistic 
tip structures, including nonsymmetric tips, is desirable for a good qualitative reproduction 
of experimental STM images. However, it is quite unlikely that the relative orientation of 
the sample surface and the local tip apex geometry in STM experiments is of high symmetry, 
which has been commonly assumed in the vast majority of STM simulations to date. Mandi 
et al. studied the effect of asymmetric relative tip-sample orientations on the STM contrast 
of the W(llO) metal surface and of the HOPG(OOOl) surface jn] employing a three- 
dimensional Wentzel-Kramers-Brillouin (3D-WKB) electron tunneling theory. It was found 
that the STM images can be substantially distorted due to tip geometry effects. A physical 
explanation was provided based on the real-space shape of the electron orbitals entering the 


orbital-dependent tunneling transmission formula in the 3D-WKB method Q , see Ea.m 
in Appendix. Motivated by the ideas of Hagelaar et al. and based on the methodology 
of Mandi et al, in the present work a new concept of obtaining information about the 
local spatial orientation of the STM tip in real instruments is introduced. The concept 
is substantiated by a combination of STM experiments and large scale simulations taking 
the HOPG(OOOl) surface. Goncomitantly, the qualitative visual analysis of STM images is 
advanced by quantifying their correspondence in terms of relative brightness correlations. 

The paper is organized as follows: The proposed correlation analysis method is introduced 
in section m followed by its application to the HOPG(OOOl) surface. We analyze and discuss 
our results in section uni and summarize our findings in section ITVl The appendix reports a 
brief summary of the 3D-WKB tunneling theory with an arbitrary tip orientation. 


II. METHOD 


To quantitatively compare the experimental (EXP) and simulated (SIM) constant-current 
pographs, the definition of the relative 

nn 

C at bias voltage I 4 is needed [UJ, : 


topographs, the definition of the relative brightness of a given two-dimensional (2D) contour 


5c(x,14) 


zc(x,Vk) - zcjy: min 5 14) 
^ci ^max5 14) - .2c( ^min? VkV 


( 1 ) 


where zc{x, I 4 ) is the apparent height of the constant-current contour C above the surface 
lateral x = Xij position at bias voltage I 4 obtained by G G (EXP, SIM}, ^^(^(xinin, 14) 
and .2c(xmax, 14 ) respectively have the smallest and largest apparent heights in the 2D scan 
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area, thus due to the definition, i?c(xmin, H;) = 0 and -Bc(xmax, Vk) = 1- Assuming that all 
Bc{xij, Vk) contours consist of x Ny points (i = 1,..., A^,, j = 1,..., Ny), the mean value 
of the relative brightness in a given bias voltage range of Ny bias values (A: = 1,..., Ny) can 
be calculated as 


_ 1 Nv Nx 


( 2 ) 


Using the same resolution of the scanning area in the experiment and in the simulations 
resulting in relative brightness contours of x Ny lateral points in both cases, it is possible 
to quantitatively compare the Sexp and Bsiu contours in the corresponding bias voltage 
range of Ny bias values by calculating their correlation coefficient as 


r = 


X 


X 


{Nv Nx 

^ ^ '^[BEXp{xij, Vk) 

k=li=lj=l 

{ Nv Nx Ny 

^ ^ '^[BEXp{xij, Vk) 

k=l i=l j=l 

{Nv Nx 

Yl ^ '^[Bsmixij, Vk) 

k=l i=l j=l 


— -SEXp][-SsiM(a^ij, 14) — -Ssim] 

I 

— -Bexp]^ 



(3) 


The (Pearson product-moment) correlation coefficient r measures the degree of linear rela¬ 
tionship between the -BEXp(a^ij, 14) and -BsiM(a^p, 14) datasets. Due to the definition, the 
values of r are bounded to the range of [-1, -|-1]. r = -|-1 corresponds to a perfect positive 
linear relationship that is desirable when comparing relative brightness contours between 
experiment and simulations. Obtaining r = -|-1 would mean that the simulation repro¬ 
duces the experimental data perfectly, r = — 1 means a perfect negative linear relationship, 
e.g., this would be the result of calculating the correlation coefficient of exactly oppositely 
corrugated contours, r = 0 means that there is no linear relationship between the contours. 

Another statistical measure for the difference between experimental and simulated con¬ 
tours is the mean squared error, 

e£i eS, E£i[Bexp(x«, 14) - U)l". 

A perfect agreement of contours is obtained at MSE=0, and it is desired that MSE is mini¬ 
mal comparing experimental and simulated contours for obtaining the best agreement. For 
selected contours and bias voltages we found good correspondence between minimal MSE 
and maximal correlation. However, MSE is not bounded from above, and this makes the 
analysis of MSE distribution and the interpretation of maximal MSE difficult. Therefore, 
we excluded using this measure in our statistical analysis. 
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The calculation of the correlation coefficient in Eq. (|2j) was presented in the more general 
case of taking 2D relative brightness contours. However, the same method can be specifi¬ 
cally applied to one-dimensional (ID) relative brightness profiles by setting Ny = 1. This 
approach will be used in the paper for the (1100) direction of the HOPG(OOOl) surface 


since experimental data m is available for such a case. To calculate the relative bright¬ 
ness correlations between the experiment and simulations, profiles shifted to start with their 
minimum value, 5c(a:i=ij=]j_l4) = 0 are taken. A detailed discussion justifying this was 
given in section 3.2. of Ref. jll|. 

Since in the simulations the tip material (TIPMAT), atomic arrangement/geometry (TIP- 
GEO), and orientation described by the Euler angles (6*o, V’o) can be chosen in practically 
infinite ways, the corresponding relative brightness profiles are dependent on these parame¬ 
ters: 

55/m(x, 14) = 55/m(x, 14, TIPMATtipgeO) ^O) t/'o), and similarly, the correlation coeffi¬ 

cient is r = r(TIPMATTiPGEO) ^0) 00) l/'o)- In fhc present work, we consider TIPMAT=W 
(tungsten) and TIPGEO G (blunt, sharp} tip models. The Wbiunt fip is represented by an 
adatom adsorbed on the hollow site of the W(llO) surface and the Wgharp Hp is modeled 
as a pyramid of three-atoms height on the W(llO) surface. More details on the used tip 
geometries can be found in Ref. [l2| . These tip models are expected to bracket the range of 
possible tip sharpnesses in experiments as extremely blunt tips with flat surfaces would pro¬ 
vide no contrast-resolution at all and sharp pyramids would likely be very unstable during 
prolonged tip scans. Moreover, carbon-contaminated tips with a C atom at the apex can be 
excluded due to a dramatic decrease of the tunneling current . 

In our simulations the three-dimensional Wentzel-Kramers-Brillouin (3D-WKB) electron 
tunneling theory with arbitrary tip orientations is employed, see Appendix, which is imple¬ 


mented in the 3D-WKB-STM code 


Q, 0, 


cessfnlly applied in a number of theoretical 


investigations 


y. 

15 


Recently, the 3D-WKB method was suc- 


-[l9j] and combined experimental-theoretical 


20 |. 


Constant-current brightness profiles are calculated along the (1100) direction (x-axis of 
Fig. [I]) containing the three characteristic positions of the HOPG(OOOl) surface: hollow (h), 
G-carbon and /3-carbon, see inset of Fig. [H The experimental averaged brightness data with 
Ny = 1 and = 46 points are taken from Fig. 4 of Ref. 0] in the interval of [-1 V, 1 V] 
with 0.1 V steps. In the simulations the current values are chosen for each corresponding 
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FIG. 1: Schematic view of the STM tip above the HOPG surface. The rotation of the 
local coordinate system of the tip with respect to that of the sample surface is described 
by the Euler angles (do, 0O) V^o)- Inset shows the positions of the characteristic h, a and (3 
sites of the HOPG(OOOl) surface along the (1100) direction. 


bias voltage in such a way that the lowest apparent height of each constant-current contour 
is 2;siM(xmin, 14) = 5.5 A (pure tunneling regime). The relative brightness profiles are 
calculated by using the introduced Wbiunt and Wgharp tip models for a set of tip orientations 
described by the Euler angles: do G [0°,30°], <po G [0°,175°], 4o G [0°,355°] with 5° steps. 
The Euler angles are visualized in Eig. [H do angle describes the rotation with respect to the 
X axis, transforming the z axis to z'. Additionally, 4o and are rotation angles around 
the z' and z axes, respectively, as Eig. [D shows. The exact meaning of the Euler angles is 


mathematical 
in Refs 


Q, 


Iv formulated in the rotation matrix in Eg. (I All 11 in Appendix and explained 
Uj. Altogether 7 x 36 x 72 = 18144 tip orientations are considered. Eor 
this selection we used the general symmetry property of the rotation matrix in Eq. flAllH : 
(do,4o)'*Ao) = (“^0)00 + 7r,4o + tt) and the mirror symmetry of the HOPG surface above 
the h — a — j3 line: (do,4o5'*Ao) = (—do, —4o, —t^o)- Gorrelation coefficients in Eq.(|3]) are 
calculated between the experimental and a large number of simulated relative brightness 
profiles in the negative (-1 V < E < 0 V, Ny = 10), positive (0 V < E < 1 V, Ny = 10) 
and full (-1 V < E < 1 V, Ny = 20) bias voltage ranges. 
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FIG. 2: Constant-current STM images illustrating the variety of observed STM contrasts 
above the HOPG(OOOl) surface in the tunneling regime for = 0o = 0°: a) hexagonal 
contrast (both a- and /3-carbons are bright; G = 1 V, = 90°), b) triangular contrast 
(only /3-carbons are bright; G = 0.1 V, = 90°), c) triangular contrast with striped 
feature {V = 0.1 V, //>o = 120°). The STM images are calculated above the shaded 
rectangular area shown in the inset of Fig. [D using the Wbiunt tip model. Inset shows the 
relative orientation of the Wbiunt tip with respect to the HOPG(OOOl) surface in each 

subfigure. 


III. RESULTS AND DISCUSSION 


We recall that the STM contrast of the HOPG(OOOl) surface can change substantially 
depending on the tunneling and tip parameters 3, y, [u, 12, 21|. A selection of the possible 
STM contrasts in the tunneling regime is shown in Fig. El Here, the two nonequivalent 
carbon atoms of HOPG (a and /3) are primarily responsible for the different STM contrasts 
[hexagonal contrast in Fig. Et) and triangular contrast in FigJ2))]. Particular rotations 
of the STM tip were shown to result in striped STM images (ll[, affecting the secondary 
contrast features [Fig. Eb)]- In the near contact regime multiple scattering effects and tip- 
sample forces also play an important role in the STM contrast appearance |22|, e.g., a shift 
of the maximum brightness from the /3-carbon to the hollow (h) position of HOPG was 


demonstrated by Ondracek et al. 


2l| . Note that we restrict our study to the pure tunneling 


regime corresponding to the used experimental data [12[ and to the validity of the 3D-WKB 
method jn]. The diversity of the observed STM contrasts above the HOPG(OOOl) surface 
surely contains information about the local geometry of the tip apex in STM measurements, 
therefore HOPG(OOOl) is an ideal candidate to illustrate the applicability of our statistical 
correlation analysis method combining large scale STM simulations with experiments. 


Fig. El shows the calculated relative brightness correlation histograms for the two consid- 
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FIG. 3: |G| < 1 V relative brightness correlation histograms calculated by using 18144 tip 
orientations for: a) Wbiunt fip, b) Wgharp tip. Part c) reports the sum of the histograms in 
a) and b). The correlation histograms for the negative, positive and full bias ranges are 
shown using Eq.(l3]) in the [0.5, 1] range with 0.001 resolution. 


ered tungsten tip models in 18144 tip orientations and the sum of the two histograms in 
Fig. |3h). The maximal correlation between the experiment and simulations is found at ap¬ 
proximately 0.97 in the negative and at approx. 0.95 in the positive bias range for both tips. 
However, we cannot conclude that the tip orientations belonging to the maximal correlation 
are the best since there is a large number of other orientations within a few percent from 
the maximum correlation well above 0.9. Analyzing the correlation distribution, it is clearly 
seen that much more tip orientations provide better correlation values in the negative com¬ 
pared to the positive bias range for both tip models. This effect is even more evident for the 
Wsharp tip, where the correlation distributions have two distinct peaks for the negative and 
positive bias at around 0.93 and 0.66, respectively. The presented statistics for the relative 
brightness correlation taking a large number of tip orientations confirm the significance of 
the findings of Ref. {ll], where the simulated brightness profiles obtained at positive bias 
for the Wsharp Hp model in high symmetry orientations resnlted in much lower correlation 
with the experiment than in the negative bias voltage range. No such large differences were 
found for the Wbiunt fip af either bias polarities. This suggests that the Wbiunt fip is more 
likely to be present in a wide range of bias voltages in the experiment than the Wsharp Hp- 


The minimal correlation between the experimental and simulated brightness profiles is 
found at 0.55 for the Wsharp tip at positive bias voltages, whereas for the Wsharp tip at 
negative bias voltages and for the Wbiunt Hp at all considered bias voltage ranges the minimal 
correlation is above 0.7. Once more, this suggests a more likely Wbiunt than Wsharp tip in 
the experiment since various local rotations of the Wbiunt tip do not give worse correlations 
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with the experiment than 0.7, whereas there are particular local rotations of the Wgharp tip 
at positive bias voltages with much worse correlations. 

The presented relative brightness correlation histograms provide information about the 
distribution of the correlation values in terms of the number of simulated tip orientations 
within a particular correlation range with the experimental brightness data. This presenta¬ 
tion of the correlation statistics, however, cannot tell which specihc tip orientations give the 
best or worst correlations with the experiment. To assign the most or least likely orientations 
of the STM tip in the experiment for the given tip model, we need another representation 
of the correlation data. Therefore, we complement our analysis by calculating correlation 
maps: r(Wbiunt, 6 * 0 ,0o, ^o) and r(Wsharp, 6 * 0 , 0o, V’o)- 

Figs, m E] and El show the calculated relative brightness correlation maps for the two con¬ 
sidered tungsten tip models in the negative, positive and full bias voltage range, respectively. 
r(0o, 'tpo) two-dimensional maps are shown as a function of Oq. Note that 9o = 0° corresponds 
to the same ^-axis of the surface and the tip, and in this case 0 o and V’o denote the same 
type of rotations around the common 2 :-axis. As a result, we obtain striped r{(/)o,'ipo) corre¬ 
lation maps for Oq = 0° [panels a) and b)]. For 0o > 0° these maps quickly change to show 
more complicated correlation distributions [panels c)-n)]. Most importantly, Figs. HI El and 
El show the most (least) likely tip orientations (dp, 0O) V^o) in the experiment in the given bias 
interval corresponding to bright (dark) regions bounded by black (white) contours within 
2% relative to the maximum (minimum) correlation value for each 9q assuming the model 
tip apex geometry. Overall, we find that the regions close to the maximal and minimal 
correlations can be differently affected by the bias range considered for the mapping for dif¬ 
ferent tip apex geometries. These results emphasize the importance of a large experimental 
dataset for reliable application of the proposed procedure. Considering the favorable and 
unfavorable orientations for the given tip models, we find that the (^OiV^o) positions of the 
indicated regions close to the maximum and minimum correlations in the r{(f)o,'ijjo) maps 
are fairly stable with respect to the change of Oq. This means that the specific ( 00 , ^o) 
Euler angles are representative for the likely (bright regions) and unlikely (dark regions) tip 
orientations in the STM experiment, irrespective of 9q. Based on our results, we find that 
the favored tip-sample relative orientations are far from being symmetric. 

We introduce the area ratios as the number of tip orientations (area) within the denoted 
regions in Figs. HJ El and El divided by the area of the r(0O)0o) niaps (36 x 72). These 
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FIG. 4: -1 V < G < 0 V negative bias range correlation analysis. Relative brightness 
correlation distributions r(do,0o,V’o) for Wbiunt Rp [first column: a), c), e), g), i), k), m)] 
and Wsharp tip [second column: b), d), f), h), j), 1), n)] for the following fixed 9o angles: 
a)-b) 0°, c)-d) 5°, e)-f) 10°, g)-h) 15°, i)-j) 20°, k)-l) 25°, m)-n) 30°. Most (least) likely tip 
orientations in the experiment in the given bias interval correspond to bright (dark) 
regions bounded by black (white) contours idithin 2% relative to the maximum (minimum) 
correlation value in each subfigure assuming the model tip apex geometry. 
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FIG. 5: 0V<G<1V positive bias range correlation analysis. Relative brightness 
correlation distributions r(do,0o,V’o) for Wbiunt tip [first column: a), c), e), g), i), k), m)] 
and Wsharp tip [second column: b), d), f), h), j), 1), n)] for the following fixed 9o angles: 
a)-b) 0°, c)-d) 5°, e)-f) 10°, g)-h) 15°, i)-j) 20°, k)-l) 25°, m)-n) 30°. Most (least) likely tip 
orientations in the experiment in the given bias interval correspond to bright (dark) 
regions bounded by black (white) contours kHthin 2% relative to the maximum (minimum) 
correlation value in each subfigure assuming the model tip apex geometry. 
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FIG. 6: |G| < 1 V full bias range correlation analysis. Relative brightness correlation 
distribntions r(do,0o,V’o) for Wbiunt Rp [first column: a), c), e), g), i), k), m)] and Wgharp 
tip [second column: b), d), f), h), j), 1), n)] for the following fixed 9o angles: a)-b) 0°, c)-d) 
5°, e)-f) 10°, g)-h) 15°, i)-j) 20°, k)-l) 25°, m)-n) 30°. Most (least) likely tip orientations in 
the experiment in the given bias interval correspond to bright (dark) regions bounded by 
black (white) contours within 2% relative tiSthe maximnm (minimum) correlation value in 
each subfignre assnming the model tip apex geometry. 

























area ratios at fixed Oq can be interpreted as the likelihood of favorable or unfavorable tip 
orientations in the experiment assuming the considered tip geometry in the given bias range. 
The area ratios alone, however, are not sufficient to identify the most or least likely tip 
orientations in the experiment since the maximum and minimum correlation values vary 
considerably depending on 6q. 

To further analyze the correlation maps in Figs. 01 ElandEl the evolutions of the maximum 
and minimum correlation values and the calculated area ratios with 9q are reported in Fig. 
[71 This figure also allows comparison between the different bias voltage ranges and the 
two considered tip models. We find that the maximum correlation is increasing and the 
minimum correlation is decreasing with increasing 6q for all bias voltage ranges. This results 
in a larger difference between the maximum and minimum correiations with increasing 6q. 
It is interesting to note that the maximum correlation values are always larger than 0.9 for 
the Wbiunt fiP) whereas this is true only in the negative bias range for the Wgharp tip. In the 
positive and full bias ranges the maximum correlation above 0.9 is achieved for 6q > 20°, i.e., 
for a much smaller number of considered tip orientations. On the other hand, the minimum 
correlation values are always smaller for the Wgharp compared to the Wbiunt tip- These 
findings cleariy suggest that the Wbiunt fip is more likely to be present in the experiment in 
an enhanced bias voltage range than the Wgharp tip. 

In Fig. [71 at negative bias voltages the two tips provide similar maximum correlation 
values as a function of 9q. In such case the area ratios can be used to decide which tip is 
more likely in the experiment since the corresponding area ratios are proportional to the 
number of tip orientations within the maximum correlation, and such larger area ratios favor 
a given tip. We find that the area ratios are generally larger for the Wbiunt compared to 
the Wsharp tip. Area ratios close to the correlation maximum mean that more orientations 
can provide better correlation values for the Wbiunt than for the Wsharp tip. On the other 
hand, area ratios close to the correlation minimum mean that more orientations provide 
correlations close to the minimum for the Wbiunt compared to the Wsharp tip. This is, 
however, not a problem in the present case since the minimum correlations are always larger 
for the Wbiunt compared to the Wsharp tip- Therefore, based on the number of favorable tip 
orientations, we can also conclude that the blunt tungsten tip is indeed more likely in the 
experiment than the sharp tip in the |V| < 1 V bias voltage range. 

In order to check the robustness of our results we performed the correlation analysis 
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FIG. 7: Analysis of the correlation maps in Fig. 0] (at negative bias), Fig. [5] (at positive 
bias) and Fig. E] (at full bias) in the |V| < 1 V bias range. Top row: The evolution of the 
maximum and minimum correlation value in the r(do, ipo, ipo) maps with 9q. Bottom row: 
The 00-evolution of the area within 2% relative to the maximum and minimum correlation 
values (respectively bounded by the black and white contours in Figs. 01 0 and Ej) in 
relation to the area of the r(0O)'*Ao) map (36 x 72). These area ratios at fixed 6q can be 
interpreted as the likelihood of favorable or unfavorable tip orientations in the experiment 
assuming the considered tip geometry. Left and right parts respectively correspond to data 

obtained by Wbiunt and Wgharp tip models. 


with simulated brightness profiles obtained by taking the contributions of four extra next- 
neighbor atoms of the tip apex atom in the tunneling current calculations using the 3D-WKB 
method. We find that the correlation maps are quantitatively very similar to those obtained 
by the one-apex tip for 0o < 20°. For larger 0o-fihmg the emergence of multiple tip apices 
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FIG. 8: |G| < 0.3 V relative brightness correlation histograms calculated by using 18144 
tip orientations for: a) Wbiunt fiP) b) Wgharp tip. Part c) reports the sum of the histograms 
in a) and b). The correlation histograms for the negative, positive and full bias ranges are 
shown using Eq.(l3]) in the [0.5, 1] range with 0.001 resolution. 

distorts the simulated brightness profiles and consequently worsens the agreement with the 
experiment, manifesting as dramatically reduced correlation values (down to 0.35 at Oq — 25° 
and 0.13 at 6q = 30°) for particular (^o^'^Ao) ranges. Based on this, we can conclude that 
our findings are robust for 6q < 20°, i.e. for a small tilting of the tip z-axis. 

To investigate the effect of the bias voltage on the obtained results, we recalculated the 
correlation statistics in the |G| < 0.3 V bias voltage range that corresponds to the low bias 
regime used in typical STM imaging experiments of HOPG. This analysis used redefined 
negative (-0.3 V < G < 0 V, Ny = 3), positive (0 V < G < 0.3 V, Ny = 3) and full (-0.3 

V < G < 0.3 V, Ny = 6) bias ranges. Fig. [S] shows the recalculated relative brightness 
correlation histograms for the two considered tungsten tip models in 18144 tip orientations 
and the sum of the two histograms in Fig. |Ht). We find qualitatively similar results as in 
the |G| < 1 V bias range reported in Fig. |21 The main differences in Fig. |H]in comparison 
to Fig. 13] are: i) there is a longer tail of the correlation distributions extending toward lower 
values for both tips, resulting in much lower minimum correlations (e.g., 0.26 for the Wgharp 
tip at positive bias voltages and 0.58 for the Wbiunt tip at all bias ranges), ii) the maximum 
correlations are increased to 0.99 at negative bias for both tips, iii) the difference between 
the two distinct peaks of the correlation distributions for the negative and positive bias in 
case of the Wgharp tip is reduced, but still significant (above 0.1). 

Fig. [9] shows the evolutions of the maximum and minimum correlation values and the 
calculated area ratios with 9q obtained from the r(do, V’o) correlation maps in the \ V\ < 0.3 

V bias voltage range. We find that the main discussed tendencies in Fig. Ejare not affected 
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FIG. 9: Extracted data from the correlation maps in the \V\ < 0.3 V bias voltage range. 

Top row: The evolution of the maximum and minimum correlation value in the 
r(do, 00,'*Ao) maps with Oq. Bottom row: The 6*o-evolution of the area ratio, for explanation 
see the caption of Fig. [71 Left and right parts respectively correspond to data obtained by 

Wbiunt and Wsharp tip models. 


in the low bias regime. However, the area ratios within 2% of the maximum correlation are 
systematically larger for the Wgharp than for the Wbiunt tip in the negative bias range. Since 
the maximum correlations are above 0.93 for for both type of tips in this bias interval, this 
suggests that more tip orientations of the Wgharp tip result in better agreement with the 
experiment than of the Wbiunt tip at low negative bias, -0.3 V < E < 0 V. The indications 
of a favored Wbiunt tip in the experiment are, however, not affected in the other considered 
low bias regimes. 

Although using larger bias ranges is better for the statistical analysis, the tip may become 
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unstable in the experiment at larger bias voltages, thus making the assignment of the tip 
geometry and orientation more difficult. In general, we suggest that the primary decision for 
the quality of the STM tip in an experiment has to be based on the comparison between the 
maximum and minimum relative brightness correlations between two (or more) tip models, 
and the secondary decisive factor should be the introduced area ratio measure that gives 
information on the number of likely or unlikely tip orientations. 


IV. CONCLUSIONS 


In scanning probe experiments the scanning tip is the source of one of the largest un¬ 
certainty as very little is known about its precise atomic structure and stability. Since the 
atomic structure and electronic properties of the tip apex can strongly affect the contrast of 
STM images, it is very difficult to experimentally obtain predictive STM images in certain 
systems. To tackle this problem we proposed a statistical correlation analysis method to ob¬ 
tain information on the local geometry and orientation of the tip used in STM experiments. 
We defined the relative brightness correlation of constant-current topographs between ex¬ 
perimental and simulated data, and analyzed it statistically for the HOPG(OOOl) surface in 
combination with two tungsten tip geometries in 18144 orientations. The simulations were 
performed using the 3D-WKB electron tunneling theory based on first principles electronic 
structure calculations. We find that a blunt tip model provides better correlation with the 
experiment for a wider range of tip orientations and bias voltages than a sharp tip model. 
A favored sharp tip is indicated at low negative bias only. From the correlation distribution 
we proposed particular tip orientations that are most likely present in the STM experiment, 
and likely excluded other orientations. Importantly, we find that the favored relative tip- 
sample orientations do not correspond to high symmetry setups that are routinely used in 
standard STM simulations. The demonstrated combination of large scale simulations with 
experiments is expected to open up the way for a more reliable interpretation of STM data in 
the view of local tip geometry effects. Moreover, the introduced correlation analysis method 
could be useful for other scanning probe imaging techniques as well. 
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APPENDIX: 3D-WKB TUNNELING THEORY 


Mandi et al. have developed an orbital-dependent electron tunneling model with arbitrary 
tip orientations for simulating scanning tunneling microscopy (STM) measurements 
within the three-dimensional (3D) Wentzel-Kramers-Brillouin (WKB) framework based on 

nnnng^ 


previous atom-superposition theories [IJ, lla, llH, l23- 


Here, this method is briefly de¬ 


scribed, which was used in the paper for the HOPG(OOOl) surface in combination with 
tungsten tips. The model assumes that electrons tunnel through a tip apex structure con¬ 
sisting of a few atoms, and transitions between individual atoms of this tip apex structure 
and a suitable number of sample surface atoins, each described by the one-dimensional (ID) 
WKB approximation, are superimposed [l3|, ll5|. Since the 3D geometry of the tunnel junc¬ 
tion is considered, the method is a 3D-WKB atom-superposition approach. The advantages, 
particularly computational efficiency, limitations, and the potential of the 3D-WKB method 
were discussed in Ref. |l^ . 

The electronic structure of the surface and the tip is included in the model by taking the 
atom-projected electron density of states (PDOS) obtained by ab initio electronic structure 
calculations j^. The orbital-decomposition of the PDOS is necessary for the description 
of the orbital-dependent electron tunneling . The energy-dependent orbital-decomposed 
PDOS functions of the ith sample surface atom with orbital symmetry a and the jth tip atom 
with orbital symmetry r are denoted by ng^(E) and respectively. In the present 

work a G {s,Py,Pz,Px} atomic orbitals for the carbon atoms on the HOPG surface and 
T G {s,Py,Pz,Px, dxy, dyz, d^z^_j. 2 , dxz, dj. 2 _y 2 } orbitals for the apex atoms of blunt and sharp 
tungsten tips are considered. The total PDOS function is the sum of the orbital-decomposed 
contributions: 


n: 


n- 


:(E) = T.n’sAE), 

a 

-(E) = y 


(A4) 

(A5) 


Note that a similar decomposition of the Green’s functions was reported within the linear 
combination of atomic orbitals (LGAO) framework in Ref. |27l |. 

Assuming elastic electron tunneling at temperature T = 0 K, the tunneling current at 
the tip apex position Rt/p and bias voltage V is given by the superposition of atomic 
contributions from the sample surface (sum over i), superposition of atomic contributions 
from the tip apex structure (sum over j) and the superposition of transitions from all atomic 
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orbital combinations between the sample and the tip (sum over a and r): 

I (R„P, V)^Y.T.T.^'Jr(RTIF,V). 


(A6) 


One particular current contribution can be calculated as an integral in an energy window 
corresponding to the bias voltage V as 

rj, (Rtip, V) = T,r {Ef, + eU, V, dp) 

X (£;| + eU) <r {e^ + el7 - eV) dU. (A7) 

Here, e is the elementary charge, h is the Planck constant, and Ep and Ep are the Fermi 
energies of the sample surface and the tip, respectively. The jh factor ensures the correct 
dimension of the electric current. The value of e has to be determined by comparing the 
simulation results with experiments, or with calculations using standard methods, e.g., the 
Bardeen approach (28|. In our simulations e = 1 eV was chosen that gives comparable 
current values with those obtained by the Bardeen method implemented in the BSKAN 


code 


m 


30(]. Note that the choice of e has no qualitative influence on the reported results. 


and a rigorous comparison between the 3D-WKB and Bardeen methods in relation to STM 
experiments performed on HOPG jl^ was reported in Ref. [nj. 

In Eq.(|MD,T.r (E,R,dp) is the orbital-dependent tunneling transmission function, and 
it gives the probability of the electron tunneling from the r orbital of the jth tip atom to the 
cr orbital of the Ah surface atom, or vice versa, depending on the sign of the bias voltage. 
The conventions of tip —>■ sample tunneling at positive bias voltage (R > 0) and sample —)■ 
tip tunneling at negative bias (V < 0) are used. The transmission probability depends on 
the energy of the electron [E), the bias voltage (V), and the relative position of the jth tip 
atom and the Ah sample surface atom (dp- = Rtip ~ Note that Rr/p corresponds 
to the position of the tip apex atom. The following form for the transmission function is 
considered [^: 

(e| + eU, V, dp-) = exp{-2n{U, V)\dij\}xUe,j, (AS) 

Here, the exponential factor corresponds to an orbital-independent transmission, where all 
electron states are considered as exponentially decaying spherical states (23|, l2a, ImI] , and it 
depends on the distance between the jth tip atom and the Ah surface atom, |dp|, and on 
the vacuum decay. 


'^s + At + gR 


-eU 


(A9) 
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For using this k an effective rectangular potential barrier in the vacuum between the sample 
and the tip is assumed, ips and ipT are the electron work functions of the sample surface 
and the tip, respectively, m is the electron’s mass and h the reduced Planck constant. The 
remaining factors of Eq. (IA8ll are responsible for the orbital dependence of the transmission. 
They modify the exponentially decaying part according to the real-space shape of the elec¬ 
tron orbitals involved in the tunneling, i.e., the angular dependence of the electron densities 
of the atomic orbitals of the surface and the tip is taken into account as the square of the real 
spherical harmonics Xcr(%,0p) and respectively. It is important to note that 

the angles are given in the respective local coordinate system of the surface (without primes) 
and the tip (denoted by primes). This distinction of the local coordinate systems is crucial 
to describe arbitrary tip orientations that correspond to a rotation of the tip coordinate 
system by the set of Euler angles (do,0o, V’o) with respect to the surface coordinate system 
ji nl |. The transformation between a vector defined in the local coordinate system of the tip, 
{x',y',z'), and a vector defined in the local coordinate system of the sample, {x,y,z), is 
given by 


fx'\ 


y 

V-V 


— ^(^0; 00; t^o) 




y 

\^J 


(AlO) 


with the rotation matrix: 


^(6>o,0o,'0o) = 

^ cos 00 cos 00 — sin 0 o sin 0 o cos Oq cos 0 o sin 0 o + sin 0 o cos 0 o cos Oq sin 0 o sin Oq ^ 


— sin 00 cos 00 — cos 0 o sin 0 o cos Oq — sin 0 o sin 0 o + cos 0 o cos 0 o cos ^o cos 0 o sin ^o 
sin 00 sin ^o “ cos 0 o sin ^o cos ^o 


(All) 


//\ //\ _ 

The polar and azimuthal angles (0T, 0)^) given in both real spherical harmonics in Eq. flASIl 

correspond to the tunneling direction, i.e., the line connecting the Ah surface atom and 
the jth tip atom, as viewed from their local coordinate systems (denoted by no prime and 
prime, respectively), and they have to be determined for each surface atom-tip atom {i — j) 
combination from the actual tip-sample geometry. A schematic view of an STM tip with 
rotated local coordinate system by the Euler angles (6*o,0o,0o) above the HOPG(OOOl) 
surface is shown in Eig. [H d, 0 and d are also shown for a given i — j pair. Eor more details 
of the 3D-WKB formalism, see Refs. Q ,Q, and for a rigorous comparison between the 
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3D-WKB and Bardeen methods in relation to STM experiments performed on HOPG, see 
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